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Abstract 

An  isothermal  three-dimensional  model  describing  mass,  momentum  and  species  transfer  in  the  cathode  of  a  proton  exchange  membrane 
fuel  cell  has  been  used  to  study  four  different  flow-distributors:  interdigitated,  coflow  and  counterflow  channels,  and  a  foam.  A  quantitative 
comparison  of  the  results  shows  that  the  interdigitated  channels  can  sustain  the  highest  current  densities,  followed  in  descending  order  by 
the  foam,  the  counterflow  and  the  coflow  channels.  The  foam  yields  the  most  uniform  current  density  distribution  at  higher  currents,  but  care 
should  be  taken  as  to  its  permeability  to  avoid  unreasonably  high-pressure  drops. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

In  view  of  ever  increasing  levels  of  environmental  pollu¬ 
tion  and  thus  a  desire  to  replace  the  fossil-fuel -based  economy 
with  a  cleaner  alternative,  the  fuel  cell  has  in  recent  years 
become  a  prime  candidate  as  a  power  source  for  transport 
and  stationary  applications.  The  potential  use  of  fuel  cells 
ranges  from  distributed  power  sources  and  portable  applica¬ 
tions,  such  as  laptops  [1]  or  even  for  the  future  dismounted 
soldier  [2],  to  vehicles. 

One  such  type  of  cell  is  the  proton  exchange  membrane 
fuel  cell  (PEMFC),  a  schematic  representation  of  which  is 
shown  in  Fig.  1 .  The  basic  cell  consists  of  two  porous  elec¬ 
trodes,  termed  the  anode  and  the  cathode,  separated  by  a 
proton  conducting  membrane.  The  porous  electrodes  are 
made  of  a  composite  material,  containing  carbon  cloth  and 
a  hydrophobic  agent,  such  as  polytetrafluorethylene.  Each 
electrode  has  a  thin  layer  containing  an  electrocatalyst,  such 
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as  platinum,  that  is  dispersed  on  the  carbon  cloth  and  is 
in  contact  with  the  membrane,  usually  a  hydrated  perfluo- 
rinated  sulfonic  acid  polymer.  In  addition  a  bipolar  plate, 
essentially  graphite  into  which  flow  and  cooling  channels 
have  been  machined,  is  situated  adjacent  to  each  electrode. 
In  the  course  of  operation,  an  oxidant,  usually  oxygen  from 
air  which  is  either  dry  or  humidified  to  some  extent,  is  fed 
at  the  inlet  on  the  cathode  side  and  transported  to  the  elec¬ 
trolyte/cathode  interface;  the  fuel  on  the  other  hand,  normally 
hydrogen,  is  fed  at  the  anode  inlet  and  is  transported  to  the 
electrolyte/anode  interface.  The  reactions  occurring  at  these 
interfaces  are  then 

2Tb -+  4H++4e“  at  the  anode,  (1) 

O2  +  4H+  — >  2H2O  at  the  cathode,  (2) 

which  are  termed  the  hydrogen  oxidation  reaction  (HOR)  and 
the  oxygen  reduction  reaction  (ORR),  respectively.  Thus,  the 
protons  produced  at  the  anode  are  transported  through  the 
membrane  to  the  cathode,  whilst  the  electrons  can  drive  a 
load  through  an  external  circuit. 

In  recent  years,  a  number  of  mathematical  models  have 
been  developed  in  an  attempt  to  understand  the  phenomena 


E.  Birgersson,  M.  Vynnycky  /  Journal  of  Power  Sources  153  (2006)  76-88 


77 


Nomenclature 

A  area  (m2) 

A  integration  constant  for  asymptotic  solution 
Aio  volumetric  exchange  current  density  (A  irf  3) 

B  transformation  tensor 

B  integration  constant  for  asymptotic  solution 

c  molar  concentration  (molm-3) 

D  diffusion  tensor  (m2  s  1 ) 

Djj  diffusion  coefficients  for  molar  diffusive  flux 
relative  to  a  molar-averaged  velocity  (m2  s-1) 
Djj  diffusion  coefficients  for  mass  diffusive  flux 
relative  to  a  mass-averaged  velocity  (m2  s-1) 
Vjj  binary  Maxwell-Stefan  diffusion  coefficients 

(m2  s-1) 

S  effective  oxygen  permeability  in  the  agglom¬ 

erates  (mol  m-1  s-1) 

A p  pressure  drop  (N  m-2) 

ex,  ey,  e-  coordinate  vectors 
E  potential  (V) 

F  Faraday’s  constant  (As  mol-1) 

ft  nucleus  effectiveness  factor  for  the  agglomer¬ 

ate  model 
h  height  (m) 

f)  relative  humidity 

i  ciurent  density  (Am-2) 

(iv)  volume  ciurent  density  (A  m-3) 

I  dimensionless  current  density 

K  permeability  tensor  (m2) 

L  length  (m) 

M  mean  molecular  mass  (kg  mol- 1 ) 

M  dimensionless  mean  molecular  mass 

Mj  molecular  mass  of  species  i  (kgmol-1) 

Mj  dimensionless  molecular  mass  of  species  i 
n  number  of  electrons  consumed  in  the  ORR  per 

oxygen  molecule 

n  unit  vector  in  the  normal  direction 

n,-  mass  flux  of  species  i  (kgm-2  s-1) 

p  pressure  (N  m-2) 

P  power  density  (W  m-2) 

r  radius  of  agglomerate  nucleus  (m) 

R  gas  constant  (J  mol- 1  K- 1 ) 

Re  Reynolds  number 

S  denominator  for  transformation  of  diffusion 

coefficients  (m2  s-1) 

t  unit  vector  in  the  tangential  direction 

T  temperature  (K) 

v  velocity  (m  s-1) 

V  volume  of  the  representative  elementary  vol¬ 

ume  (m3) 
w  width  (m) 

wi  mass  fraction  of  species  i 


W  transformation  tensor 

Xi  molar  fraction  of  species  i 

X,Y  dimensionless  coordinates 

X  transformation  tensor 

Greek 

a  coefficient  for  water  transport  by  electro¬ 

osmosis  in  the  membrane 
ar  cathodic  transfer  coefficient  for  the  ORR 

y  porosity 

S,j  Kronecker  delta 

A  =  l /(Rea2)  dimensionless  parameter  for  asymptotic 
solution 

A  dimensionless  fraction  of  species  i  for  asymp¬ 

totic  solution 
ij  overpotential  (V) 

k  permeability  (m2) 

pt  dynamic  viscosity  (kg  m- 1  s- 1 ) 

§  stoichiometry 

p  density  (kg  m-3) 

a  =  hf/L  dimensionless  number 

<tsU|  standard  deviation  for  current  density  (A  m-2) 

T  parameter  for  current  density  expression 

4>  general  tensor 

0  —  (2  +  4q,).Mh20  —  Afo,  dimensionless  number 
ft?  =  At  a  dimensionless  number 

Subscripts 

0  equilibrium  (reference) 

a  active  layer 

cell  cell 

p  porous  backing 

pol  polymer  electrolyte  in  the  active  layer 

f  flow  channel 

FFO  water 

02  oxygen 

N2  nitrogen 

avg  average 

Superscripts 

cat  catalytic  region 

g  gas 

in  inlet 

out  outlet 

ref  reference 

vap  vaporization 

Miscellaneous  symbols 
()  superficial  average 

(}®  intrinsic  average 

[]  scale 
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Fig.  1 .  Schematic  of  a  proton  exchange  membrane  fuel  cell. 


occurring  in  a  PEM  fuel  cell.  Since  a  complete  fuel  cell  model 
would  have  to  address  mass,  momentum,  species  and  heat 
transfer  in  gas  and  liquid  phases  in  a  three-dimensional  geom¬ 
etry,  as  well  as  the  electrokinetics  for  the  ORR  and  HOR, 
most  models  choose  to  focus  on  only  some  of  these  aspects 
at  a  time.  The  first  models  to  appear  were  one-dimensional 
followed  by  two-dimensional;  see  [3]  for  a  list  of  these 
models.  Lately  three-dimensional  models,  based  on  compu¬ 
tational  fluid  dynamics  (CFD),  have  also  begun  to  appear 
[4-18].  The  main  differences  between  the  3D  models  are 
the  flow-distributors  they  study,  whether  or  not  multiphase 
flow  is  modelled  and  whether  or  not  a  uniform  tempera¬ 
ture  is  assumed  throughout  the  cell.  Most  of  these  three- 
dimensional  CFD  models  consider  one-phase  flow  conditions 
[4-6,8-12,15,17,18],  i.e.  only  take  the  gas-phase  into  account 
in  the  flow  fields  and  porous  backings.  While  this  is  a  valid 
approximation,  provided  the  cell  temperature  is  sufficiently 
high  and  the  partial  pressure  of  the  water  vapor  remains  below 
the  saturation  value,  water  vapor  tends  to  condense  in  the  cell 
and  form  a  liquid  phase,  especially  at  higher  current  densities. 
Berning  et  al.  [4]  and  Zhou  and  Liu  [5]  developed  models  to 
predict  the  local  behaviour  of  species,  temperature  and  cur¬ 
rent  density  in  a  PEM  fuel  cell  equipped  with  parallel  flow 
channels;  Berning  et  al.  [4]  later  used  theirs  for  a  parameter 
study  of  operational,  geometric  and  material  properties  [6]. 
An  extension  by  Berning  and  Djilali  [7]  to  account  for  the 
liquid  water  phase  allowed  for  predictions  of  the  liquid  sat¬ 
uration  distribution  in  the  cell.  Dutta  et  al.  [8]  examined  the 


performance  of  the  PEM  fuel  cell  for  straight  channels  with 
and  without  the  porous  backings  as  well  as  the  water  transport 
through  the  membrane.  The  model  was  later  applied  by  Dutta 
et  al.  [9]  to  elucidate  the  impact  of  serpentine  flow  channels 
on  cell  performance.  The  model  by  Jen  et  al.  [10]  was  used 
to  predict  cell  performance  with  straight  channels.  Kumar 
and  Reddy  [11]  examined  the  impact  of  geometrical  features 
such  as  channel  width,  depth  and  shape  for  serpentine  flow 
fields  on  cell  performance  and  later  also  porous  metal  foams 
[12],  for  which  they  observed  a  more  uniform  current  den¬ 
sity  distribution  compared  to  a  multi-parallel  channel  flow 
field  design.  Lee  et  al.  [13]  extended  the  model  developed 
by  Shimpalee  and  Dutta  [14]  to  study  serpentine  flow  fields, 
including  experimental  validation  for  the  overall  water  bal¬ 
ance.  Mazumder  and  Cole  derived  both  a  gas-phase  [15]  and 
a  multi-phase  model  [16]  to  predict  cell  performance  and 
amount  of  liquid  water.  Senn  and  Poulikakos  [17]  found  that 
a  porous  foam  led  to  improved  mass  transfer  and  higher  cell 
performance  compared  to  straight  and  serpentine  flow  chan¬ 
nels.  Shimpalee  and  Dutta  [14]  developed  a  model  to  study 
various  parameters,  such  as  species  and  temperature  distri¬ 
butions  for  straight  channels.  The  interdigitated  flow  channel 
design  was  compared  with  straight  channels  by  Um  and  Wang 
[18],  who  found  that  the  interdigitated  can  sustain  a  higher 
mass-transport  limited  current  density  than  the  straight  chan¬ 
nels.  More  complex  flow  field  distributors  have  also  been 
studied,  such  as  the  alternative  flow  field  design  comprising 
tree  network  channels  [19].  Based  on  a  ID  model,  Senn  and 
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Poulikakos  [19  ]  found  that  the  tree  network  structure  can  sus¬ 
tain  higher  currents  than  a  parallel  flow  field  design,  provided 
that  the  fuel  cell  is  designed  in  a  shape  that  is  reminiscent  of 
a  double-staircase. 

To  extend  the  work  on  three-dimensional  modelling,  we 
have  conducted  a  quantitative  comparison  of  the  performance 
of  four  common  flow-distributors:  parallel  channels,  run  in 
both  coflow  and  counterflow,  interdigitated  channels  and  a 
porous  distributor,  such  as  a  foam;  comparison  is  carried  out 
with  respect  both  to  cell  performance  and  current  density 
uniformity  over  the  catalytic  layer.  To  limit  the  computa¬ 
tional  requirements,  we  assume  that  the  anode  side  is  run 
at  such  conditions  that  it  is  able  to  fully  sustain  any  current 
created  at  the  cathode,  i.e.  a  fully  moisturized  anode  that  is 
run  at  high  stoichiometry.  In  addition,  we  assume  that  suffi¬ 
cient  cooling  is  provided  to  keep  the  cathode  isothermal,  a 
not  unreasonable  preliminary  assumption  in  view  of  the  tem¬ 
perature  distributions  obtained  by  [4,14],  who  show  that  the 
temperature  only  varies  by  a  couple  of  degrees  in  the  cathode. 
These  assumptions  enable  us  then  to  consider  isothermal, 
three-component,  gas-phase,  three-dimensional  laminar  flow 
in  the  flow-distributor  and  adjacent  porous  backing  on  the 
cathode  side  only. 

The  mathematical  model,  consisting  of  mass,  momentum 
and  species  transport  equations,  as  well  as  the  geometries  of 
the  flow-distributors  considered,  are  introduced  in  Section  2. 
We  focus  also  on  how  a  rather  detailed  agglomerate  model 
for  the  electrochemical  aspects  of  the  active  layer,  derived 
by  [20],  can  be  simply  implemented  into  the  present  formu¬ 
lation.  Details  of  the  numerical  solver  used,  CFX-4.4  [21], 
followed  by  its  verification  against  an  asymptotic  solution 
obtained  previously  [3]  are  given  in  Section  3.  The  results 
from  different  flow-distributors  are  then  compared  and  dis¬ 
cussed  in  Section  4.  We  finish  with  conclusions  in  Section  5. 


2.  Model  description 

2.1.  Flow-distributor  geometry 

The  electrochemical  reactions  that  occur  at  the  active  lay¬ 
ers  depend  on  a  sufficiently  fast  transport  of  reactants  to,  and 
products  away  from,  the  active  sites  so  as  to  limit  concen¬ 
tration  overpotentials.  Towards  this  end,  the  bipolar  plates 
contain  grooved  channels,  which  can  take  a  number  of  dif¬ 
ferent  shapes.  Amongst  the  most  common  designs  today  are: 

(a)  parallel  channels,  with  only  one  pass  over  the  porous 
backing,  run  in  coflow,  as  shown  in  Fig.  2a; 

(b)  parallel  channels,  with  only  one  pass  over  the  porous 
backing,  run  in  counterflow,  as  shown  in  Fig.  2a; 

(c)  interdigitated  channels,  where  channels  are  terminated, 
in  order  to  force  the  flow  into  the  porous  backing,  see 
Fig.  2b; 

(d)  a  porous  material,  such  as  a  foam;  here,  the  entire  surface 
of  the  porous  backing  is  in  contact  with  the  gas  flow,  see 


Fig.  2.  Schematic  of  the  flow-distributors  considered:  (a)  parallel  channels, 
which  can  be  run  in  coflow  or  counterflow;  (b)  interdigitated  flow  channels; 
(c)  foam. 


Fig.  2c,  in  contrast  to  the  channel  based  flow-distributors, 
which  contain  ‘dead’  zones  between  the  channels.  Note 
that  we  consider  inlets  and  outlets  of  the  same  size  as 
for  (a-c),  even  though  a  flow  field  comprising  a  porous 
medium  can  be  operated  with  continuous  inlets  and  out¬ 
lets.  The  difference  in  pressure  drop  and  cell  performance 
between  the  inlet  and  outlet  configuration  chosen  here  and 
a  continuous  counterpart  is  negligible  due  to  the  length 
to  width  aspect  ratio  L/w>f^>  1  (see  Table  1); 

(e)  serpentine  flow  channels,  comprising  one  long  channel 
with  many  passes  over  the  porous  backing; 

(f)  a  combination  of  some  of  the  above,  e.g.  (a  and  e), 
(b  and  e). 

We  focus  here  on  (a-d). 

2.2.  Governing  equations 

2.2.1.  Flow  channels 

In  the  flow  channels,  we  solve  for  the  momentum  and 

continuity  of  mass,  given  by 

V  •  (pv)  =  0,  (3) 


V  ■  (pv  <S>  v) 


=  -V 


2 

3 


/iV'V 


+V  •  (/z((Vv)  +  (Vv)T)), 


(4) 


where  v  is  the  velocity,  p  the  density,  p  the  pressure  and  /i  is 
the  dynamic  viscosity.  The  transport  equations  for  the  ternary 
gas  mixture,  comprising  oxygen,  water  and  nitrogen,  are 


V  • 


wo2 

Wh20 


=  V  • 


V»o2 

Vwh20 


(5) 
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Table  1 


Base-case  parameters 


Physical  parameters 

Yb 

0.3 

Yf 

0.99 

®02,H20 

3.98  x  10~5  m2  s"1  at  363  K,  1  atm  [20] 

®02,N2 

2.95  x  10”5  m2  s"1  at  363  K,  1  atm  [20] 

®h2o,n2 

4.16  x  10-5  m2  s"1  at  363  K,  1  atm  [20] 

Eo 

1.18V 

Mo , 

3.2  x  10-2  kgmoP1 

Mh^O 

1.8  x  10_2kgmol_1 

M N, 

2.8  x  10_2kgmol_1 

10-12m2 

Kf 

10“ 10  m2 

n 

1.812  x  10-5  kgm-1  s— 1 

F 

96487  As  mol-1 

a 

0.3 

R 

8.3 14  J  mol-1  K"1 

m 

10-2  kg  moP1 

[pi 

1 .0  kg  m-3 

vap 

Ph2o 

4.7  x  104  N  m-2  [27] 

Operating  conditions 

pout 

1  atm 

L 

0.1  m 

30% 

T 

80  °C  (353  K) 

hf 

10-3  m 

h. 

3  x  10~4  m 

w\ 

5  x  10-4  m 

w 

2  x  10“3  m 

Agglomerate  model  parameters 

r 

5  x  10“7  m 

ar 

0.78 

Aia 

3  x  103  Am-3 

» 

10-11  molnP1  s_1 

K 

10“5  m 

Ka 

0.3 

Tpol 

0.3 

n 

4 

cref 

c02 

6.2  mol  m-3  at  fj  =  30%,  1  atm  and  353  K 

1  Wf  corresponds  to  the  half-width  of  the  flow  channels  or  inlet  for  the  foam 
in  the  computational  domains  due  to  symmetry  conditions. 


where  w o2  and  ii>h,o  are  the  mass  fractions  of  oxygen  and 
water  and  D  is  the  diffusion  tensor. 

2.2.2.  Porous  backing/foam 

For  porous  regions,  we  have  to  define  superficial  and 
intrinsic  properties.  Superficial  averages  are  defined  as 

(4«>  =  ^,f  <MV,  (6) 

v  Jv 

and  intrinsic,  as 

<«<E>  s  v®  X4,v-  ,7) 

where  V  is  the  total  volume  of  the  representative  elemen¬ 
tary  volume  (REV)  and  ]/g)  is  the  volume  of  the  gas  in 
the  REV.  With  the  porosity,  y  =  V*g)/V  the  two  averages  are 


related  through 

(<t>)  =  /(^>(g).  (8) 

Conservation  of  mass  and  momentum  is  given,  respectively, 
by 

V  •  ((p>(g)  (v»  =  0,  (9) 


V-«p>(g)  (v)  <g>  (v»  +  /xK_1  ■  (v> 


-l 


=  -V  (  (P)ig)  +  ■  (V] 

3  y 


+V  ■  (j(V  (v)  +  (V  (v»T) 

where  K  is  the  permeability  tensor. 

The  species  transport  equations  are  described  by 

«)) 


(10) 


V-  f<p}(g)<v)  ( 

V  V(wH2o) 

=  V  •  f(p)(sV(D>(g) 


V<iuo2)(g) 

V(wh2o)(S) 


(11) 


where  (D)®  is  the  total  mass  diffusion  tensor,  containing 
contributions  from  an  intrinsic  effective  mass  diffusion  tensor 
and  an  intrinsic  hydrodynamic  dispersion  tensor.  For  a  more 
detailed  discussion  of  these,  see  [3], 


2.3.  Boundary  conditions 


2.3.1.  Inlet 

At  the  inlet  (see  Fig.  3),  we  prescribe  the  normal  velocity 
and  the  gas  composition  for  the  channel  distributors: 

v  •  ev  =  f/m,  v  •  ey  =  v  ■  ez  =  0,  wo2  =  w o2> 

wh2o  =  Wh2o;  (12) 

in  addition,  for  the  counterflow-distributor,  we  require 

v  ■  e.Y  =  -t/in,  v  •  ey  =  v  ■  e-  =  0,  wo2  =  i"o,. 

wH20  =  wH-)0’  (12) 

for  the  second  channel  (see  Fig.  3b); 

For  the  foam: 

(v)  ■  eA  =  C/in,  (v)  •  ey  =  (v)  ■  e-  =  0, 

<wo2)(g)  =  wiq2,  <wh2o)(S)  =  wS2o-  (I4) 


2.3.2.  Outlet 

At  the  outlet,  we  specify  the  pressure  and  the  stream- 
wise  gradients  of  the  velocities  and  species  are  set  to  zero, 
corresponding  to  fully  developed  flow  conditions.  For  the 
channels: 


(15) 
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Fig.  3.  Schematic  of  the  computational  unit  cells:  (a)  parallel  channels  run  in  coflow;  (b)  parallel  channels  run  in  counterflow;  (c)  interdigitated  channels;  (d) 
foam. 


(ev  •  V)(v  •  ev)  =  0, 

(16) 

(e.v  •  V)(v  •  e-)  =  0, 

(17) 

e  t  -  V  w07  =  eA  •  V  ujh2o  =  0; 

(18) 

correspondingly,  for  the  foam: 

(P>(g)  =  P°at, 

(19) 

(e.v  •  V)((v>  •  ev)  =  0, 

(20) 

(e,  •  V)«v>  •  e,)  =  0, 

(21) 

e*  •  V<wo2)(g)  =  eA  •  V<utH2o)(g)  =  0. 

(22) 

2.3.3.  Walls 

At  the  walls  of  the  channels,  we  specify  no-slip,  no  normal 
flow  and  no  componental  flux: 

v  =  0,  (23) 

n  •  Vu)o2  =  n  •  VwhiO  =  0,  (24) 

where  n  is  the  unit  normal  to  a  wall.  In  the  porous  backing  and 
foam,  zero  shear  stress,  no  normal  flow  and  no  componental 
flux  conditions  are  applied;  these  are,  respectively 


n  •  (v)  =  0,  (29) 

n  •  V(wo2)(g)  =  n  ■  V(wH2o)<g)  =  0;  (30) 

correspondingly,  for  the  channels: 

(n  •  V)(v  •  t)  =  0,  (31) 

n  •  v  =  0,  (32) 

n  •  V wq2  =  n  •  V u)h,o  =  0.  (33) 


2.3.5.  Channel/porous  backing  interface 

At  the  interface  between  the  porous  backing  and  the  flow 
channels,  we  couple  the  pointwise  velocities,  normal  and 
shear  stresses  in  the  plain  fluid  (flow  channels)  with  their 
superficial  counterparts  in  the  porous  medium.  The  mass  frac¬ 
tions  and  fluxes  of  oxygen  and  water  are  continuous  across 
the  interface. 

2.3.6.  Active  region/porous  backing  interface 

The  active  region  in  the  cathode  is  sufficiently  thin  to  allow 
us  to  treat  it  as  a  boundary  condition  for  the  porous  backing. 
Using  Faraday’s  Law,  the  mass  fluxes  of  oxygen  and  water 
can  be  found  as 


(n  •  V)((v>  ■  t)  =  0, 
n  ■  (v)  =  0, 


(25) 

(26) 


(no2)  •  ey  = 


Mq2  (i)  •  ev 
4  F 


(34) 


n  •  V<W02)(g)  =  n  •  <VWH2o)(g>  =  0,  (27) 

where  t  is  the  unit  tangent  to  a  wall. 

2.3.4.  Symmetry  conditions 

The  flow-distributors  considered  in  this  paper  are  all  of 
periodic  character,  see  Fig.  2,  allowing  us  to  reduce  the  com¬ 
putational  domain  for  each  by  introducing  unit  cells,  see 
Fig.  3,  with  the  appropriate  symmetry  conditions.  For  the 
foam  and  porous  backing: 

(28) 


(nH2o)  • ev  = 


(1  +  2qQMh2o  (i)  •  ey 
2 F 


(35) 


where  (no2)  and  (nn2o)  are  the  mass  fluxes,  a  the  amount 
of  water  dragged  through  the  membrane  with  each  proton, 
Mj  the  molecular  mass  of  species  i,  (i)  the  superficial  current 
density  and  F  is  Faraday’s  constant.  The  superficial  velocity 
from  the  reaction  can  be  derived  from  the  theory  of  multi- 
component  mass  transfer  as 


(p>(g)  (v)  •  ev 


•  e 


4  F 


(n  •  V)((v>  ■  t)  =  0, 


—((2  +  4a)A/n2o  —  Mq2). 


(36) 
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2.4.  Constitutive  relations 


Assuming  an  ideal  gas,  the  gas  density  can  be  found  from 


pM 

~RT’ 


(37) 


where  M  —  (wo2/M02  +  wh2o/Mh2o  +  wn2/MN2)_1  is 
the  mean  molecular  mass,  R  the  universal  gas  constant  and  T 
is  the  temperature.  The  mass  fraction  of  nitrogen  is  given  by 


wn2  =  1  —  wo2  —  »h2o- 


(38) 


We  note  here  that  the  molar  fractions  x ;  are  related  to  the 
mass  fractions  Wj  by 


X)  = 


WjM 

Mt  ' 


(39) 


Furthermore,  in  the  most  general  case,  the  dynamic  viscos¬ 
ity,  pt,  is  also  a  function  of  the  composition,  but  for  simplicity 
it  is  considered  constant  in  this  paper.  The  corresponding  con¬ 
stitutive  properties  in  the  porous  backing  and  foam  are  the 
same,  but  based  on  intrinsic  values. 

An  expression  for  the  diffusion  tensor  [22]  can  be  found 
from  the  diffusion  coefficients  Dp  for  the  molar  diffusion 
flux,  relative  to  a  molar-averaged  velocity  frame,  as 


On  = 

77o2,n2(xo27?h2o,n2  +  (1  —  -*o2)X>o2,h2o) 

(40) 

5 

D\2  — 

xo27?h2o,n2(7?o2,n2  —  T>o2,i[2o) 

(41) 

5 

Ih\  — 

xh2o77o2,n2(77h2o,n2  —  7702,h2o) 

(42) 

5 

The  cross  terms  in  the  diffusion  tensor  D  are  around  one  to 
two  orders  of  magnitude  lower  for  the  operating  parameters 
in  this  study  than  the  diagonal  terms,  allowing  us  to  neglect 
their  contributions,  whence 


0 

D22 


(50) 


In  the  porous  media,  i.e.  in  the  foam  and  porous  backing,  we 
apply  a  Bruggeman  correlation  for  the  superficial  effective 
mass  diffusion  tensor: 


(D)  =  y3/2  D. 


(51) 


The  permeabilities  of  the  porous  backing  and  foam  are  taken 
to  be  isotropic 

K  =  KSij,  (52) 

where  k  is  the  permeability. 


2.5.  Electrokinetics  and  the  active  layer 

An  expression  is  still  required  for  the  current  density  (i)-ey 
given  in  equation  (36).  A  novel  feature  of  this  paper  is  to 
implement  an  expression  from  an  agglomerate  model  derived 
by  Jaouen  et  al.  [20],  and  subsequently  validated  experimen¬ 
tally  [24].  Although  [20]  considers  an  active  layer  of  finite 
thickness,  it  is  possible  to  demonstrate  that,  for  our  purposes, 
the  active  layer  need  not  be  resolved  explicitly,  but  rather  can 
be  treated  implicitly  as  a  boundary  condition.  Some  further 
details  are  as  follows. 

The  volumetric  current  density  (iv),  given  by  [20],  is 
approximated  as 


72hiO,n2(-«h2o77oo,No  +  (1  -  *&o)£,02,h2o) 

Oil  =  - - - - - ,  (43) 

S  —  xo2£>h2o,n2  +  xh2o77o2,n2  +  xn2T>o2,h2Oi  (44) 


{iv)  =  Ar'0(l  -  TpolXl  -  Ka)  exp 


3 


K>w 


c 


ref 

O, 


(53) 


where  'Dp  are  the  binary  Maxwell-Stefan  diffusion  coef¬ 
ficients.  Since  we  use  the  mass  diffusion  flux  relative  to 
the  mass-averaged  velocity,  the  following  transformation  is 
required: 


D  =  BWX-1DXW-1B-1 , 

(45) 

„  „  (,  » >nXj\ 

B  —  $ij  U^i  I  1  )  5 

V  XnWjJ 

(46) 

X  —  XiSp,  i,  j  —  1,2, 

(47) 

3 

II 

S 

0-3 

II 

to 

(48) 

where  Sy  is  the  Kronecker  delta.  Here,  i,  j  =  1  correspond  to 
O2,  i,  j  =  2  correspond  to  H2O,  and  n  =  3  corresponds  to  N2. 
The  binary  Maxwell-Stefan  diffusion  coefficients  are  cor¬ 
rected  for  pressure  and  temperature  via  [23]: 

/  rji  \  3/2 

VijiT,  P)=j[y0)  Vij(T°’ Po1  m 


where  Aio  is  the  volumetric  exchange  current  density  in 
the  agglomerates,  ypo|  the  volume  fraction  of  the  poly¬ 
mer  electrolyte  in  the  agglomerate  nucleus,  (co,/c)  = 

(ioo,)®  {p)^ / Mo2  the  molar  concentration  of  oxygen,  ar 
the  cathodic  transfer  coefficient  for  the  ORR,  t)  the  overpo¬ 
tential  at  the  cathode  (defined  negative)  and  y.d  is  the  volume 
fraction  of  pores  in  the  active  layer.  $  is  the  nucleus  effec¬ 
tiveness  factor,  defined  as 

3  (  1  1 

Tr  \tanh(7>)  Tr 

with  T  given  by 

Aipjl  -  yPoi)  exp 

nFD  ’  C  ’ 

where  D  is  an  effective  oxygen  permeability  in  the  agglomer¬ 
ates,  n  the  number  of  electrons  consumed  in  the  ORR  per  oxy¬ 
gen  molecule  and  r  is  the  radius  of  the  agglomerate  nucleus. 
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This  agglomerate  model  was  validated  by  [24]  for  a  small 
PEM  fuel  cell  with  an  area  of  2  cm2  at  conditions  well  above 
the  stoichiometric  flow  rate,  allowing  the  cell  there  to  be  mod¬ 
elled  one -dimensionally,  since  concentration  gradients  in  the 
streamwise  direction  in  the  cathode  could  be  neglected.  In 
this  paper,  however,  the  aim  is  to  capture  and  study  two-  and 
three-dimensional  effects.  The  one-dimensional  agglomerate 
model  can  be  used  for  this  purpose  by  noting  that  the  geome¬ 
try  of  the  active  layer  studied  is  slender,  i.e.  its  thickness,  h-d, 
is  much  smaller  than  its  width  (w)  and  breadth  (L),  implying 
that  only  gradients  in  the  normal  direction  in  the  active  layer 
will  contribute  to  mass  transfer;  consequently,  this  allows  us 
to  treat  the  active  layer  locally  as  one-dimensional.  Although 
we  omit  the  details  here,  it  is  possible  to  show  by  scale  analy¬ 
sis  that  we  will  not  have  any  gradients  of  oxygen  in  the  normal 
direction  in  the  active  layer,  since  the  magnitude  of  the  dif¬ 
fusion  coefficient  in  the  agglomerates  is  (?(10_13)m2  s  , 
as  compared  to  <9(10-5)  m2  s-1  in  the  pores.  The  additional 
approximation  of  reducing  the  active  layer  to  a  boundary 
condition  corresponds  to  an  infinite  effective  proton  conduc¬ 
tivity  in  the  active  layer.  The  total  current  density  is  then  given 
locally  by  (i)  ■  ev  =  (i)  =  (iv)  ha. 

Jaouen  et  al.  [20]  discerned  four  different  regimes,  where 
the  Tafel  slope  doubles  or  even  quadruples,  and  subsequently 
supplied  the  experimental  validation  to  support  these  [24]. 
In  regime  1,  the  active  layer  is  controlled  by  Tafel  kinet¬ 
ics  and  is  first  order  in  the  oxygen  concentration.  Regime 
2  displays  a  doubling  of  the  Tafel  slope,  due  to  the  active 
layer  being  governed  by  Tafel  kinetics  and  oxygen  diffu¬ 
sion  in  the  agglomerates,  but  still  remains  first-order  in 
the  oxygen  concentration.  A  doubling  of  the  Tafel  slope 
is  observed  in  the  third  regime,  where  the  active  layer 
is  controlled  by  the  Tafel  kinetics,  in  addition  to  proton 
migration.  The  oxygen  dependence  here  is  half-order.  The 
final  regime,  the  fourth,  shows  a  quadrupling  of  the  Tafel 
slope,  and  is  attributable  to  an  active  layer  controlled  by 
Tafel  kinetics,  proton  migration  and  oxygen  diffusion  in  the 
agglomerates.  The  oxygen  dependence  is  half-order,  as  in 
regime  3. 

By  assuming  that  we  have  no  resistance  to  proton  migra¬ 
tion,  we  limit  the  validity  of  the  expression  to  regimes  1 
and  2,  i.e.  to  an  active  layer  controlled  by  Tafel  kinet¬ 
ics  at  low  overpotentials,  and  at  higher  overpotentials  by 
Tafel  kinetics  coupled  with  oxygen  diffusion  resistance  in 
the  agglomerates.  We  are  thus  able  to  capture  the  dou¬ 
bling  of  the  Tafel  slope  due  to  mass  transfer  limitations 
in  the  agglomerates,  although  the  doubling  due  to  proton 
migration  resistance  and  a  quadrupling  of  the  slope  are  not 
modelled. 

Finally,  we  note  that  of  more  use,  as  regards  judging  fuel 
cell  performance,  than  the  overpotential, is  the  cell  voltage, 
£ceii,  where 

1]  =  £cell  -  Eq,  (56) 

and  Eq  is  the  equilibrium  potential. 


3.  Numerics  and  verification 

A  commercial  computational  fluid  dynamics  code,  CFX- 
4.4,  based  on  finite  volumes,  was  used  to  implement  the 
model  outlined  above.  As  can  be  seen  in  Fig.  2,  all  of  the  flow- 
distributors  are  periodic  in  the  spanwise  direction.  Hence,  a 
representative  computational  unit  cell,  with  symmetry  bound¬ 
ary  conditions  on  both  sides  in  the  spanwise  coordinate,  can 
be  introduced.  The  unit  cells  for  the  distributors  are  chosen  so 
that  the  area  of  the  active  layer  and  the  contact  area  between 
the  channels  and  porous  backing  are  the  same  for  all,  with 
the  exception  of  the  foam,  which  covers  the  entire  porous 
backing  as  depicted  in  Fig.  3. 

Using  a  structured  mesh,  each  flow-distributor  with  porous 
backing  was  resolved  as  follows: 

(a)  Parallel  channels  run  in  coflow :  The  channel  contained 
104  computational  cells  and  the  porous  backing  2  x  104 
cells,  giving  a  total  of  3  x  104  cells.  Mesh  indepen¬ 
dence  was  assured  by  comparing  with  a  mesh  comprising 
1.2  x  105  cells  and  the  difference  was  found  to  he  ~  1  % 
for  the  average  current  density. 

(b)  Parallel  channels  run  in  counterflow'.  The  two  channels 
with  104  cells  each  and  the  porous  backing  with  4  x  104 
cells,  giving  a  total  of  6  x  104  cells.  Mesh  indepen¬ 
dence  was  assured  by  comparing  with  a  mesh  comprising 
2  x  1 05  cells  and  the  difference  was  found  to  be  ~  1  %  for 
the  average  current  density. 

(c)  Interdigitated  channels'.  The  two  channels  with  104  cells 
each  and  the  porous  backing  with  4  x  104,  giving  a  total 
of  6  x  104  cells.  Mesh  independence  was  assured  by 
comparing  with  a  mesh  comprising  2  x  105  cells  and  the 
difference  was  found  to  be  ~4%  for  the  average  current 
density.  At  higher  currents,  i.e.  ( i )  >  4  A  cm-2,  more  than 
6  x  104  cells  were  necessary  to  resolve  the  flow  in  the 
porous  backing.  We  did  not  pursue  these  higher  current 
densities,  however. 

(d)  Foam:  The  foam  and  porous  backing,  each  with  2  x  104 
cells,  giving  a  total  of  4  x  104  cells.  Mesh  indepen¬ 
dence  was  assured  by  comparing  with  a  mesh  comprising 
1.6  x  105  cells  and  the  difference  was  found  to  he  ~  1  % 
for  the  average  current  density. 

The  computations,  carried  out  on  a  500  MHz  Compaq 
AlphaServer  with  3  GB  RAM,  required  about  1-2  h  for  lower 
current  densities  and  about  12  h  for  higher. 

The  numerical  code  was  verified  via  asymptotic  solutions 
[3],  who  showed  that  for  a  slender  two-dimensional  geometry, 
consisting  of  a  flow  channel  adjacent  to  the  porous  backing, 
closed  form  solutions  could  be  found  for  A  =  I /(Pea2)  1, 
where  Re  (=pUm  U/i)  is  the  Reynolds  number,  and  o  =  hf/L . 
with  hf  and  L  as  the  height  and  length  of  the  flow-distributor. 
The  geometry  is  depicted  in  Fig.  4  in  dimensional  form.  For 
the  asymptotic  analysis,  the  governing  equations  are  scaled 
so  that  X = x/L.  Y  =  y/hf.  For  a  current  density  which  can  be 
written  non-dimensionally  on  the  form  I  —  ,  where 
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Fig.  4.  The  two-dimensional  cathode  with  a  flow-distributor  above  the 
porous  backing. 


1=  {/)/[/]  and  [;]  is  the  current  density  scale,  the  leading  order 
solution  was  found  to  be  given  by 


Alo  /fo2(X)(<2Xo2(0)+l)\ 
OgUo2(0)(^o2(X)+i)j 

-(A0  -  B)  ( - fo2(X)  -  fo2(0) - \ 

’\mo2(X)+l)mO2(0)+l  )J 

C2X 

=  4po(0)(^?o2(0)+l)’ 

where 


A  = 


1 


Mn2 

(2 


(A4h2o  -  A4n2) 

I-  4a)£o2(0)  +  Ch2o(0) 
^o2(0)  +  l 


-  1 


(57) 


Fig.  5.  Verification  of  the  CFX-4.4  code.  ( — )  Corresponds  to  analytical 
solutions  and  markers  are  for  the  CFX-4.4  solutions;  (+)  water  mass  fraction; 
(x)  oxygen  mass  fraction. 

2D  calculation  in  CFX-4.4  were  2000.  The  mass  fractions 
obtained  with  CFX-4.4  agree  well  with  the  ones  obtained 
analytically. 


4.  Results  and  discussion 


B  = 


1 


(Mq2  -  A4n2)  +  (A4h2o  -  A4n2) 


Mn2 

'  <£Xh2o(0)  ~  (2  +  4a) 
<2X0X0)+  1 


(58) 


For  each  flow-distributor,  simulations  were  carried  out  for 
six  different  cell  potentials,  Ecen.  However,  rather  than  spec¬ 
ifying  the  inlet  velocity,  Um,  we  specify  instead,  as  tends 
to  be  the  practice  in  experimental  work  with  fuel  cells,  the 
stoichiometry,  £,  which  is  defined  by 


and 


Co2  = 


*o2 

M  ’ 


(59) 


/A.n»g2dA  _  pUinw^Ain 

!a ca,  (n$)  dA  ~~  Mo2/4+/ACat  (i)  d A'’ 


(62) 


£h2o 

=  (1  +  <ZXq2(X)Kh2o(0)  +  (2  +  4oQ(fo2(0)  -  fo2(X)) 

4Xo(0)+l 

(60) 

with,  A ij  —  Mj/[M],  M.  =  M/[M\ ,  where  [M]  is  a  molec¬ 
ular  weight  scale,  Q  =  and  0  =  (2  +  4a)Mu^o  ~ 

A4ot  ■  For  the  agglomerate  model  applied  here,  we  obtain  the 
current  density  scale: 

,  (  arF  \  /ra[p] 

M  =  Ai0(l  -  ypoi)(  1  -  Ka) exp  -  —  ri  f  , 

V  RT  )  cg;[M] 

(61) 

whence  I  =  (xo2)<d>  as  required  for  the  solution  above. 

Comparison,  shown  in  Fig.  5,  was  carried  out  in  terms 
of  mass  fraction  profiles  along  the  cathode.  All  parameters 
from  the  base  case  were  used,  see  Table  1,  but  with  A  =  960, 
<j=  ICC3  and  £2=  28.2.  The  number  of  cells  used  for  this 


where  n q2  and  ylol/  are  t^le  mass  fluxes  of  oxygen  into 
the  cathode  and  of  the  oxygen  being  consumed  at  the  active 
layer,  respectively,  and  Acat  and  Ain  are  the  total  areas  of 
the  active  layer  and  the  inlet.  This  formulation  implies  that 
the  inlet  velocity  is  iterated  for,  its  value  depending  on  the 
current  density  obtained  at  the  catalytic  layer  on  the  previous 
iteration.  We  also  specify  the  relative  humidity  at  the  inlet, 
5)ln,  given  by 


f)ln  = 


vin  „in 

ah2o  P 

vap  » 

Fh2o 


(63) 


where  p^_ PQ  is  the  vapor  pressure  of  water;  furthermore,  since 
x02An2  =  21/79  and  xq2  +  xj^o  +  xn2  =  1,  we  have  the 
inlet  compositions  for  a  given  relative  humidity  as 


xS2o  =  X 


vap 

in  P H20 


+)2 


1  —  rin 
1  xH2Q 

1  +  79/21 ' 


(64) 
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Fig.  6.  Polarization  curves  for  the  different  flow-distributors  at  stoichiometry 
§=1.5,  3,  5:  coflow  channels  (OX  counterflow  channels  (V),  interdigitated 
channels  (x),  foam  (+). 

Consequently,  x^q  and  Xq,  must  also  be  updated  at  every 
iteration. 

Earlier  work  [25,26],  albeit  in  2D  models,  suggests  that  the 
interdigitated  design  should  be  capable  of  yielding  higher  cur¬ 
rent  densities  than  either  one  of  the  straight  channel  designs 
considered  here.  Also,  we  expect  that  either  the  foam  or  the 
interdigitated  channels  give  rise  to  the  highest  average  cur¬ 
rent  densities,  since  the  interdigitated  design  forces  the  flow 
into  the  porous  backing,  whereas  the  foam  covers  the  whole 
surface  of  the  porous  backing.  The  foam  has  already  been 
shown  to  be  more  efficient  in  terms  of  higher  current  den¬ 
sities  and  power  density  than  serpentine  flow  channels  [17]. 
Parallel  channels  run  in  counterflow  might  be  expected  to  per¬ 
form  better  than  channels  in  coflow,  as  the  former  allows  for 
alternating  inlets  and  outlets,  thus  providing  for  an  increased 
mass  transfer  in  the  spanwise  direction.  This  is  indeed  the 
case,  as  can  be  seen  in  Fig.  6,  where  the  polarization  curves 
for  the  stoichiometries  1.5,  3  and  5  are  shown.  The  interdigi¬ 
tated  design  allows  for  average  current  densities  of  4  A  cm-2 
for  cell  potentials  between  0.3  and  0.4  V,  depending  on  stoi¬ 
chiometry.  Current  densities  ranging  from  2  to  3.2  A  cm-2  are 
obtainable  with  the  foam.  The  counterflow  gives  somewhat 
higher  current  densities  than  the  coflow  design,  especially  for 
the  higher  stoichiometry. 

Fig.  7  depicts  the  power  density  for  the  flow-distributors. 
In  the  lower  ranges  of  current  density,  0-0.4  A  cm-2,  the 
power  density  is  independent  of  the  type  of  distributor  used. 
As  the  current  density  is  increased,  differences  between  the 
flow  designs  become  readily  apparent,  with  the  interdigitated 
sustaining  the  highest  power  densities,  followed  by  the  foam. 
The  counterflow  channels  perform  marginally  better  than  the 
coflow  channels. 

The  performance  of  a  fuel  cell  is  judged  not  only  on  the 
magnitude  of  current  density  that  can  be  generated,  but  also 
on  the  uniformity  of  the  current  distribution  at  the  active  layer. 


Fig.  7.  Power  density,  P,  for  the  different  flow-distributors  at  stoichiometry 
§=  1.5,  3,  5:  coflow  channels  (OX  counterflow  channels  (V),  interdigitated 
channels  (x),  foam  (+). 

since  uniformity  is  linked  to  catalyst  utilization  and  degrada¬ 
tion.  As  a  measure  of  the  uniformity,  we  compare  the  standard 
deviation  of  the  current  density  for  each  flow  design,  erstd, 
defined  by 

where 

<65) 

as  illustrated  in  Fig.  8.  For  all  flow-distributors,  the  distribu¬ 
tion  becomes  more  uniform  as  the  stoichiometry  is  increased. 
At  a  stoichiometry  of  1 ,  all  the  oxygen  that  enters  the  cathode 
would  be  consumed;  by  increasing  the  stoichiometry,  the  flow 


Fig.  8.  Standard  deviation  of  the  current  density,  astd,  for  the  different  flow- 
distributors  at  stoichiometry  §=1.5,  3,  5:  coflow  channels  (OX  counterflow 
channels  (V),  interdigitated  channels  (x),  foam  (+). 


Fig.  9.  The  local  current  density  distribution  for  the  coflow-distributor  at  the 
active  layer  for  different  potentials  at  stoichiometry  £  =  3:  (a)  Eceii  =  0.782  V; 
(b)  Eceii  =  0.682  V;  (c)  Eceii  =  0.582  V;  (d)  Eceii  =  0.482  V;  (e)  Ece u  =  0.382  V; 
(f)£ceii  =  0.282  V. 


Fig.  10.  The  local  current  density  distribution  for  the  counterflow-distributor 
at  the  active  layer  for  different  potentials  at  stoichiometry  £  =  3:  (a) 
£cen  =  0.782  V;  (b)  Eceii  =  0.682  V;  (c)  Eceii  =  0.582  V;  (d)  Ece ;U  =  0.482  V;  (e) 
Eceii  =  0.382  V;  (f)  £ceii  =  0.282  V. 


becomes  less  depleted  of  oxygen,  and  is  hence  able  to  sustain 
higher  current  and  a  more  uniform  distribution.  The  paral¬ 
lel  channels,  run  in  co-  and  counterflow,  exhibit  the  highest 
deviations  at  current  densities  above  MJ.4  A  cm-2,  with  the 
coflow  being  the  less  uniform  of  the  two.  At  current  densities 
below  1.5-2  A  cm-2,  depending  on  stoichiometry,  greatest 
uniformity  is  obtained  for  the  interdigitated  design,  but  the 
extent  of  the  non-uniformity  increases  with  the  current;  ulti¬ 
mately,  the  foam  design  gives  the  greatest  uniformity.  These 
findings  are  also  reflected  in  Figs.  9-12,  where  the  local  cur¬ 
rent  density  for  different  cell  potentials  is  given.  All  distribu¬ 
tors  exhibit  a  more  non-uniform  current  density  as  the  over¬ 
potential  is  increased,  i.e.  as  the  cell  voltage  decreases.  For 
coflow  (see  Fig.  9),  the  areas  under  the  flow  channels  can  sus¬ 
tain  local  higher  current  densities  than  those  under  the  “rib”  of 
the  bipolar  plate,  where  mass  transfer  becomes  increasingly 
limiting  with  increasing  overpotential.  The  cathode  operated 
in  counterflow  behaves  similarly  (Fig.  10),  but  delivers  a 
somewhat  higher  overall  performance,  as  the  counterflow 
arrangement  allows  for  exiting  air  with  a  lower  oxygen  con¬ 
centration  in  one  channel  to  be  balanced  by  incoming  fresh  air 
in  the  two  adjacent  channels.  For  the  foam  (Fig.  11),  the  local 
current  density  is  a  function  of  the  streamwise  coordinate 
only,  although  there  are  minor  inlet  and  exit  effects,  which 
become  somewhat  more  pronounced  at  higher  current  densi¬ 
ties.  This  simplification  is  attributable  to  the  inherent  charac¬ 
teristic  of  the  foam  to  cover  the  surface  of  the  porous  backing, 
in  contrast  to  grooved  channels  in  a  bipolar  plate  comprising 
alternating  regions  of  channels  and  ribs.  Finally,  the  interdigi¬ 
tated  arrangement  (Fig.  12)  displays  an  increasingly  spanwise 


behaviour  for  the  local  current  density  at  higher  overpoten¬ 
tials.  The  overall  flow  increases  for  a  given  stoichiometry 
at  higher  current  densities,  leading  to  increased  forced  con¬ 
vective  flow  in  the  porous  backing,  whence  the  local  current 
density  becomes  more  even  in  the  spanwise  direction. 


xio5 


xlO3 


(e) 
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Fig.  1 1 .  The  local  current  density  distribution  for  the  foam  distributor  at  the 
active  layer  for  different  potentials  at  stoichiometry  £  =  3:  (a)  Eceii  =  0.782  V; 
(b)  Ecdi  =  0.682  V;  (c)  Eceii  =  0.582  V;  (d)  Eceii  =  0.482  V;  (e)  Eceii  =  0.382  V; 

(f)  ECeii  =  0.282  V. 
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Fig.  12.  The  local  current  density  distribution  for  the  interdigitated  flow- 
distributor  at  the  active  layer  for  different  potentials  at  stoichiometry  §  =  3: 
(a)  £ceii  =  0.782  V;  (b)  £ceii  =  0.682  V;  (c)  Ece M  =  0.582  V;  (d)  £cen  =  0.482  V. 


Keeping  the  pressure  drop,  A p  (=  pm  —  pout),  at  a  mini¬ 
mum  is  of  interest  in  terms  of  reducing  operating  costs  for  the 
fuel  cell,  whence  a  good  flow-distributor  should  be  able  to 
sustain  high  even  current  densities,  whilst  keeping  the  pres¬ 
sure  drop  to  a  minimum.  The  foam,  as  shown  in  Fig.  13, 
requires  the  highest  pressure  drop  to  drive  the  flow;  this  can  be 
attributed  to  the  rather  low  permeability  chosen  for  the  foam 
in  this  study.  An  increase  in  permeability  to  10-8  m2  would 
lead  to  a  reduction  of  the  pressure  drop  by  approximately  two 
orders  of  magnitude,  as  can  be  estimated  from  Darcy’s  law. 
The  pressure  drop  for  the  interdigitated  distributor  is  higher 
than  for  the  coflow  and  counterflow,  which  is  to  be  expected 
since  the  flow  is  being  forcibly  driven  through  the  porous 
backing,  in  this  study  with  a  permeability  of  10“ 12  m2.  The 
lowest  pressure  drops  are  obtained  with  the  parallel  chan¬ 
nels,  with  no  discernible  difference  in  the  magnitude  of  the 
pressure  drop  between  the  two. 


Fig.  13.  Pressure  drop,  A p,  for  the  different  flow-distributors  at  stoichiom¬ 
etry  §=1.5,  3,  5:  coflow  channels  (OX  counterflow  channels  (V),  interdigi¬ 
tated  channels  (x),  foam  (+). 


Fig.  14.  Obtained  inlet  velocity,  Uin,  for  the  different  flow-distributors  at 
stoichiometry  §=1.5,  3,  5:  coflow  channels  (OX  counterflow  channels  (V), 
interdigitated  channels  (x),  foam  (+). 


Since  we  do  not  specify  the  inlet  velocity,  but  rather  iterate 
on  it  to  obtain  a  given  stoichiometry  for  a  specified  cell  poten¬ 
tial,  it  is  of  interest  to  see  how  the  inlet  velocity  changes  with 
current  density  and  stoichiometry;  this  is  shown  in  Fig.  14. 
The  inlet  velocity  is  almost  proportional  to  the  average  cur¬ 
rent  density,  and  only  starts  to  deviate  from  linear  when  the 
density  at  the  inlet  changes  due  to  the  pressure  drop  obtained. 
The  reason  why  the  inlet  velocities  for  the  interdigitated  chan¬ 
nels  are  higher  is  because  the  inlet  area  is  half  as  large  as  for 
the  other  flow-distributors. 


5.  Conclusions 

A  study  of  different  flow-distributors,  based  on  a  gas-phase 
model  for  mass,  momentum  and  species  transport  in  the  cath¬ 
ode  of  a  PEM  fuel  cell,  has  been  considered.  Special  attention 
was  given  to  the  treatment  of  the  active  layer,  for  which  an 
agglomerate  model  developed  previously  [20]  was  used.  The 
numerical  code  used  was  verified  with  an  analytical  solution, 
also  developed  previously  [3].  The  novel  features  of  this  study 
are:  ( 1 )  a  direct  quantitative  comparison  of  the  performance  of 
different  flow-distributors  and  (2)  the  use  of  a  new  agglomer¬ 
ate  model  which  agrees  well  with  experimental  observations 
for  a  small  fuel  cell  [24] .  The  complete  validation  of  the  gas- 
phase  model  considered  here  would  require  more  detailed 
information  about  the  local  current  density  distribution;  in 
future,  this  will  be  obtainable  via  experiments  with  a  seg¬ 
mented  cell. 

The  aim  of  the  study  was  to  compare  the  performance  of 
different  flow-distributors  for  a  given  cell  at  a  given  poten¬ 
tial,  in  terms  of  four  different  quantities:  the  obtained  average 
current  density,  power  density,  standard  deviation  of  the  cur¬ 
rent  density  distribution  and  pressure  drop.  The  results  show 
that  the  interdigitated  flow-distributor  can  sustain  the  high- 
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est  current  densities,  but  at  a  higher  pressure  drop  than  the 
counterflow  and  coflow  channels.  Furthermore,  to  function 
properly,  the  interdigitated  channels  would  have  to  be  in 
contact  with  the  porous  backing  in  such  a  way  that  channel¬ 
ing  effects  are  kept  at  a  minimum;  given  the  high  velocities 
required,  even  the  slightest  gap  might  lead  to  most  of  the  flow 
going  through  the  gap  and  not  through  the  porous  backing, 
with  a  resulting  loss  of  power  density.  A  foam  distributor 
is  able  to  give  the  lowest  standard  deviation  for  the  current 
at  high  current  densities,  but  care  should  be  taken  as  to  its 
permeability  to  avoid  an  unreasonably  high  pressure  drop. 

The  present  work  was  limited  to  gas -phase  flow  and 
isothermal  conditions.  Future  work  will  seek  to  carry  out  a 
comparative  study  that  incorporates  both  the  possible  produc¬ 
tion  of  liquid  water  at  the  catalytic  layer  and  non-isothermal 
effects. 
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